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Abstract 

The site-diluted Ising model has been investigated using an improved micro- 
canonical algorithm from Creutz Cellular Automaton. For a microcanonical al- 
gorithm, the basic problem is to estimate the correct temperatures using average 



' values of the kinetic energy in the simulations of site-diluted Ising model. In this 

o ■ 

, study, the average kinetic energy has been re-described with an expression depen- 



dent on dilution x = 1 — p. The values of the temperature have been calculated 
using the new expression and the critical temperatures have been estimated from 
the peaks of specific heat for each value of dilution x. The obtained phase transi- 
tion line {kTc{p) / J , p) is in good agreement with functional prediction for the 
site-diluted Ising model. The simulations were carried out on a square lattice with 
periodic boundary conditions. 

Keywords: site-diluted Ising model;critical behavior; critical temperature; cel- 
lular automaton; square lattice; microcanonical 

1 . Introduction 

In the recent years, the Ising model has been applied with success in many 
different physical situations such as the site-diluted ferromagnet [1 — 7], the 
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pure ferromagnet[8 — 10], microemulsion[ll — 13], structural and magnetic 
phase transition[14 — 19]. Especially, the influence of defects (site or bond 
dilution) on spin systems has been reviewed by many theoretical, numerical 
and experimental investigations. The systems with defect are modelled using 
modified versions of the most popular models for pure systems, such as Ising 
and Potts models. The site-dilution on the physical systems has a significant 
effect on the critical behavior of the phase transitions. Therefore, the site- 
diluted Ising model has been applied to the problems as the determination 
of the critical temperatures for the order-disorder phase transitions [20 — 
22] . While two-dimensional pure Ising model was solved exactly many years 
ago, the diluted Ising model has not been solved yet for two and upper 
dimensions. However, it has been investigated using different simulation and 
approximation techniques such as Monte Carlo [1 — 5, 21], series expansion [7], 
renormalization group [23, 24], transfer matrix [6] and mean field [25]. The 
Hamiltonian for the site-diluted Ising model may be written as 



where S'j = ±1 is the spin at site i and the J parameter is the ferromagnetic 
spin-spin interaction energy constant. The sum is carried out over all nearest- 
neighboring (nn) spin pairs on a two-dimensional square lattice. The £j is 
called the occupation coefficients of site i. The value of £j is 1, if the spin is 
present or if the spin is absent at site i. The values of occupation coefficients 
( Ei) are randomly distributed on the lattice. Their configurational average 
takes a value in the interval of < p = (Sj) < 1 . On the other hand, the 




(1) 



<ij> 
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dilution of system is obtained hy x — 1 — p. The dilution x is for the pure 
Ising model. The aim of this work is to study the critical behavior of the 
site-diluted Ising model. In particular, we would like to locate the phase 
transition line Tc{p) and to study the thermodynamic quantities at around 
this line. The functional prediction for phase transition line Tc{p) has been 
given by 

Tc(p) = [tanh-^ (e-i-45(p-Pc))]-i (3) 

where Pc — 0.593 [3, 26]. Recently, the phase transition line of the site-diluted 
two-dimensional Ising model has been obtained using different Monte Carlo 
algorithms [1 — 5]. 

In the present paper, we have studied two-dimensional site-diluted Ising 
model using a cooling algorithm improved from Creutz cellular automa- 
ton(CCA), and wc have obtained the {kTc{p)/J, p) phase diagram. The 
CCA algorithm is a microcanonical algorithm interpolating between the con- 
ventional Monte Carlo and molecular dynamics techniques on a cellular au- 
tomaton, and it was first introduced by Creutz for the pure spin- 1/2 Ising 
model[27]. The improved CA algorithms for the versions of the Ising model in 
two and higher dimensions have been proven to be successful in producing the 
values of the static critical exponents and the critical temperatures [28 — 34]. 
In these algorithms, the temperature is not an input parameter and its value 
is obtained from the average kinetic energy of system. The estimation of 
the correct temperature for the site-diluted system using the microcanonical 
algorithms is an important problem due to dilution. The dilutions leads to 
the formation of new states for kinetic energy. Therefore, the average kinetic 
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energy expression should be rewritten depending on dilution ( a: ). In this 
study, we first investigated the dependence on dilution of the average kinetic 
energy. Furthermore, to estimate the phase transition line Tc{p), we com- 
puted the temperature variations of the specific heat {C/k) and the internal 
energy ([/) for different x — 1 — p values on the L x L square lattice with 
periodic boundary conditions. 

The paper is organized as follows. The details of the model are given in 
Section 2, the expression of average kinetic energy is modified for site-diluted 
systems in section 3, the data are analyzed and the results are discussed in 
section 4 and a conclusion is given in Section 5. 

2 . Model 

In the improved CA algorithm for site-dilute Ising model, the four vari- 
ables are associated with each site of the lattice. The value of variables on 
the each site is determined from its value and those of its nearest- neighbors 
at the previous time step. The updating rule, which defines a deterministic 
cellular automaton, is as follows: the first two variables on each site, are the 
Ising spin Bi and the occupation coefficient Si. Their values may be or 1. 
The Ising spin energy for the site-diluted model is given by Eq.l. In Eq.l 
Si — Bi — 1 — ^1. The third variable is for the momentum variable conjugate 
to the spin ( the demon ). The kinetic energy associated with the demon, 
Hfc, is an integer, which equals to the change in the Ising spin energy for 
any spin flip and its value lies in the interval (0, m). The upper limit of the 
interval, m, is equal to 16J. The total energy 

H^Hi + Hk (3) 
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is conserved. 

The fourth variable provides a checkerboard style updating, and so it al- 
lows the simulation of the Ising model on a cellular automaton. The black 
sites of the checkerboard are updated and then their color is changed into 
white; white sites are changed into black without being updated. The up- 
dating rules for the spin and the momentum variables are as follows: For a 
site to be updated, its spin is flipped and the change in the Ising spin energy 
,dHi, is calculated. If this change in energy is transferable to or from the 
momentum variable associated with this site, such that the total energy H 
is conserved, then this change is done and the momentum is appropriately 
changed. Otherwise, the spin and the momentum are not changed. 

For a given total energy, the temperature of the pure Ising model is ob- 
tained from the average value of kinetic energy, which is given by [27]: 



where n is equal to the possible kinetic energy value of a site on the spin 
system. The values of n can be zero and multiples of four for the pure Ising 
model. The expectation value in Eq. 4 is an average over the lattice and the 
number of time steps. Because of the third variable, the algorithm requires 
two time steps to give every spin of the lattice a chance to change. Thus, 
in comparison to ordinary Monte Carlo simulations, two steps correspond to 
one full sweep over the system variables. 

The cooling algorithm is divided into two basic parts, the initialization 
procedure and the taking of measurements [2 8]. In the initialization proce- 
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dure, firstly, all the spins in the lattice sites take the ferromagnetic ordered 
structure (tt) and the occupation coefficients ( Ei) are randomly distributed 
on the lattice corresponding to value of dilution x. The kinetic energy per 
site is equal to the maximum change in the Ising spin energy for the any spin 
flip using the third variables. This configuration is run during the 10.000 
cellular automaton time steps. In the next step, the last configuration in 
the disordered structure has chosen as a starting configuration for the cool- 
ing run. Rather than resetting the starting configuration at each energy, it 
is convenient to use the final configuration at a given energy as the starting 
point for the next. During the cooling cycle, energy has been subtracted from 
the spin system through the third variables {Hk) after the 200.000 cellular 
automaton steps. 

3. The modification of the average kinetic energy expression for 
site-dilute systems 

As is known, the temperature of system is obtained using average kinetic 
energy for a microcanonical algorithms. For the pure Ising model, the values 
of the kinetic energy (n) can be zero and multiples of four. But, its values 
can be zero or multiples of two dependent on x for site-diluted Ising model. 
Therefore, to estimate the temperature of system from the average value of 
the kinetic energy is a basic problem on the site-diluted Ising model simula- 
tions with a microcanonical algorithm. At the same time, the possible values 
of n at each sites of the lattice vary depending on the values of dilution x. 
If the least one of the nearest neighbor sites of a site is empty, the values of 
n can be zero or multiples of two for this site. Otherwise, n takes the values 
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zero or multiples of four. The dilution causes the new kinetic energy states 
with different frequency in the site-diluted system compared to the the pure 
Ising model. Therefore, the average value of the kinetic energy can not be 
defined accurately using Eq.(4) and the system temperature can not be esti- 
mated for the site-dilute Ising model. As a result of these, the average kinetic 
energy expression should be a function of dilution x. As is well known, the 
kinetic energy in the pure Ising model obeys the Boltzmann's law for any 
given temperature. The probability P{n), that Hk has the energy n, is 

p-nJ/kT 

Pin) = (5) 

For the site-diluted Ising model, the expressions of the partition function 
should have a different multiplier {p{x)) dependent on dilution x for each 
kinetic energy states as given below, 

Z-j;p,(x)e-^^A^. (6) 

z=0 

In this case, the expressions of the probability and the average kinetic energy 
can be re- written as a function of x. 

F{n) = eMl^ (7) 

where 2; = 0, 1, 2, I and I is equal to 8 . 

Depending on the change in the Ising spin energy for a spin flip, the 
values of kinetic energy can be equal to the values of n = 0, 4, 8, 12, 16 
for the pure system(x — 0) and the values of n = 0, 2, 14, 16 for the 



site-diluted system {x — 0.5) at any site. Here, the 2, 6, 10 and 14 values of 
kinetic energy occurs because of site-dilutions. For any value of x, the p^(x) 
can be expressed as 

p,{x) ^ {1 - xy-' x' (9) 

where k = mod{z, 2), p^ix) is equal to 1 for even z and for odd z at x — > 
for the pure Ising model. However, its values is equal to {1 — x) for even z 
and X for odd z in the site-diluted Ising model. 
3. Results and discussion 

The thermodynamic quantities of site-diluted Ising model are calculated 
using cooling algorithm. The values of the quantities are averages over the 
lattice and over the number of time steps (200.000) with discard of the first 
20.000 time steps during which the cellular automaton develops. The sim- 
ulations have been performed 10 times for initial configurations with the 
different total energy at each values oi x — 1 — p — 0, 0.05, 0.10, 0.15, 0.20, 
0.25 and 0.30 on the finite lattices with hnear dimensions L — 60, 80, 100 
and 120. 

To determine the correct temperature value corresponding to the average 
kinetic energy of system, we firstly calculated the probabilities P{n)* for 
each n value from simulations on the finite size lattice with linear dimension 
L = 120. The probabihty of each n value for kinetic energy is calculated by 

P{ny = NjNccAS, (10) 

where Nn is the number of appearance for n value on kinetic energy during 
the simulation and Nccas is the total number of Creutz cellular automaton 
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steps. In order to determine the appropriate expression for the average kinetic 
energy, P{n)* has been compared with the P{n) and F{n) probabihties in 
the Eq.(5) and the Eq.(7). The variation of the probabilities are ilhistrated 
for three different temperature values { T < Tc, T c^ii Tc and T > Tc ) 
in Fig.l. It can be seen from Fig.l that the simulation results are in good 
agreement with values of the F{n) probability in Eq.(7) for T < Tc, T ^Tq 
and T > Tq ■ At the same time, the calculated values of the probabilities 
are also in good agreement with F{n) in the interval < x < 0.4 for the 
all {Hk) values. This case shows that the temperature of the site-diluted 
system can be estimated with the expression dependent on x in Eq.(8) for 
average kinetic energy. 

The temperature dependence of the internal energy ([/) and the specific 
heat {C/k) at various values of dilution x are shown in Fig. (2) and Fig. (3) for 
single realizations of the site-dilution on a finite lattice with linear dimension 
L — 120. The values of the specific heat {C/k) and the internal energy {U) 
calculated by taking the average over the time step and lattice using the 
following formulas. 

U={-JY, eiSjSiSj )/2L' (11) 

<i3> 

C/k = L^{{U^)-{Uf)/{kTf (12) 

In these figures, the values of temperatures are estimated using Eq.(8) for 
the average value of the kinetic energy. Also, the simulations have been 
performed 10 times with various initial configurations for the each values 

9 



of the linear dimension L and of the dilution x. As seen in Fig. (2), Ising 
energy displays a continuous behavior for all dilutions and the functional 
behaviors arise in region of low temperature with increasing dilution. In 
Fig. (3), the specific heat C/k exhibits a strong peak around the critical 
temperature ior x — 0. The peak height decreases rapidly while the dilution 
increases. The functional behavior of specific heat corresponds to second- 
order phase transition. At x > 0.10 region, the specific heat shows a shoulder 
at temperatures above Tc as seen in Monte Carlo simulations [4, 5]. 

The finite-size lattice critical temperatures are also obtained from the 
specific heat maxima T^{L,x) and the infinite lattice critical temperatures 
T^'{oo,x) are estimated by analyzing the data within the framework of the 
finite size scaling theory. The critical temperature values are average of the 
estimated critical temperature values for the 10 different simulations. An 
overall error of about 3% is estimated for the values of critical temperatures. 
According to the finite size scaling theory [32], the infinite lattice critical 
temperature Tc{oo,x) is given by, 

re(oo, x) = Tc{L, x) + aL-^/", (13) 

where = 1 is the exponent associated with the divergence of correlation 
length in infinite lattice for x = 0[30]. The values of the finite lattice critical 
temperature against 1/L are illustrated for various values of dilution x in 
Fig.4. For the all dilution values, the straight lines which fit these data give 
infinite lattice critical temperature as 1/L — > 0. The critical temperature 
values for the finite lattices in Fig. 4 (a), (b) and (c) are estimated using 
n = 0,2, ...,16 in Eq.(8), n - 0,2, ...,16 and n = 0,4, ...,16 in Eq.(4) for the 
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average value of the kinetic energy, respectively. The estimated infinite lattice 
critical temperatures obtained using Eq.(4) and Eq.(8) and the functional 
prediction for phase transition line Tc{p) arc shown in Fig. 5. It can be 
seen in Fig. 5 that the functional prediction is in good agreement with the 
estimated infinite lattice critical temperatures obtained using the expression 
of the average kinetic energy dependent on x for the phase transition line 
Tc{p = 1 — a;). However, the critical temperature values estimated using n 
values for pure Ising system are quite compatible in the interval < a; < 0.05 
with the predicted phase transition line. But this compatibility disappears 
for the increased values of dilution ( x — > xc )• On the other hand, the 
estimated critical temperatures using n = 0, 2, 14, 16 values in Eq.(4) 
close to the predicted transition line only for high values ( x — > xc ) of 
dilution. The simulations results showed that one can calculate the true 
temperature of diluted system using Eq.(8) for average kinetic energy. 
4. Conclusion 

The site-dilute Ising model has been simulated using the CA cooling 
algorithm for the square lattice with the ferromagnetic interactions. The 
simulations show that the dilution leads to the formation of the new states 
for kinetic energy in microcanonical algorithm of Ising model. Therefore, the 
expression for the expected value of the kinetic energy should have a different 
multiplier {p{x)) depending on dilution x for each kinetic energy state on the 
site-diluted Ising model. The correct temperature values have been obtained 
by the expression depend on dilution (x). The infinite lattice critical tem- 
peratures T^{oo,x) are estimated by analyzing the finite-size lattice critical 
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temperatures obtained from the specific heat maxima within the framework 
of the finite size scahng theory. For phase transition fine Tc{ p = I — x), the 
estimated infinite lattice critical temperatures using expression dependent on 
X of the average kinetic energy are in good agreement with the functional 
prediction. As a result of calculations of probability -P(n), the suggested ex- 
pression (Eq.(8)) has been quite successful on the estimation of phase tran- 
sition line for site-dilute Ising model. In the simulations of the site-diluted 
Ising model with a microcanonical algorithm, the system temperature can 
only be estimated with the expression dependent on x for the average kinetic 
energy. This expression can be used for a accurate temperature account in 
simulation by a microcanonical algorithm of a system containing the defect. 
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Figure Captions: 

Figure 1. The n dependence of the probability P{n) at (a)T > Tc (b) 
T ~ Tc, (c) T < Tc for L = 100. 

Figure 2. The temperature variation of internal energy {Hi) at the various 
X values for L = 120. 

Figure 3. The temperature variation of specific heat {C/k) at the various 
X values for L = 120. 

Figure 4. The plots of the finite lattice critical Tc{L,x) temperature 
against 1/L at various x values, (a) for n = 0, 2, 14, 16 in Eq.(8), (b) for 
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n = 0, 2, 14, 16 in Eq.(4) and (c) for n = 0, 4, 8, 12, 16 in Eq.(4) 

Figure 5. The phase diagram of the site diluted Ising model on square 
lattice. The infinite lattice critical temperatures Tc{oo,x) obtained for n = 
0, 2, 14, 16 in Eq.(8), n = 0, 2, 14, 16 in Eq.(4) and n = 0, 4, 8, 12, 16 in 
Eq.(4). 
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